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1. Introduction. 

Fluid-Structure Interactions (FSI) are of great relevance in many fields of applied sci- 
entific and engineering disciplines. A comprehensive study of such problems remains 
a challenge and justifies the attention made over the last decades to propose efficient 
and robust numerical methods. We refer to [24] where different numerical procedures 
to solve the FSI problems are reviewed. One classification of FSI solution procedures 
can be based upon the treatment of the meshes with conforming or non-conforming 
mesh methods. For the first ones, meshes are conformed to the interface where the 
physical boundary conditions are imposed [26, 40, 42]. As the geometry of the fluid 
domain changes through the time, re-meshing is needed, what is excessively time- 
consuming, in particular for complex systems. 

In the present paper, we are interested in non-conforming mesh methods with a fic- 
titious domain approach where the mesh is cut by the boundary. Most of the non- 
conforming mesh methods are based upon the framework of the immersed methods 
where force-equivalent terms are added to the fluid equations in order to represent 
the fluid-structure interaction [34, 29]. Many related numerical methods have been 
developed, in particular the popular distributed Lagrange multiplier method, intro- 
duced for rigid bodies moving in an incompressible flow [19]. In this method, the fluid 
domain is extended to cover the rigid domain where the fluid velocity is required to 
be equal to the rigid body velocity. This constraint is enforced by using distributed 
Lagrange multipliers, which should be approximated on a mesh covering the structure 
and sufficiently coarse with respect to the mesh used for the fluid velocity, in order to 
satisfy the inf-sup condition. 

More recently extended Finite Element Method (XFEM) introduced by Moes, Dol- 
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bow and Bclytschko in [31] (see a review of such methods in [15]) has been adapted 
to FSI problems in [16, 30, 47, 10]. The idea is similar to the fictitious domain / 
Lagrange multiplier method mentioned above, but the fluid velocity is no longer ex- 
tended inside the structure domain and its equality with the structure velocity is 
enforced by a Lagrange multiplier only on the fluid-structure interface. One thus gets 
rid of unnecessary fluid unknowns and moreover one easily recovers the normal trace 
of the Cauchy stress tensor on the interface. We note that XFEM has been originally 
developed for problems in structural mechanics mostly in the context of cracked do- 
mains, see for example [22, 32, 44, 48, 46]. The specificity of the method is that it 
combines a level-set representation of the geometry of the crack with an enrichment 
of a finite clement space by singular and discontinuous functions. Several strategies 
can be considered in order to improve the original XFEM. Some of these strategies 
are mathematically analyzed in [25, 9]. 

In the context of fluid-structure interactions, the difficulty that present the appli- 
cations of such techniques lies in the choice of the Lagrange multiplier space used 
in order to take into account the interface, which is not trivial due to the fact that 
the interface cuts the mesh [3]. Indeed, the multiplier space, besides having good 
approximation properties, should satisfy an uniform inf-sup condition (similarly to 
more traditional fictitious domain methods [18]). In a straightforward discretization, 
it implies that the mesh for the multiplier should be sufficiently coarse in comparison 
with the mesh for the primal variables. Thus, the natural mesh given by the points 
of intersection of the interface with the global mesh cannot be used directly. An al- 
gorithm to construct a multiplier space satisfying the inf-sup condition is developed 
in [3], but its implementation can be difficult in practice. It may be thus preferable 
to work on the natural, easily constructible mesh, as outlined above. This is achieved 
in a stabilized version of the method proposed in [22] (an extension to the contact 
problems in elastostatics is also available in [23]). In the present paper, we are inter- 
ested in extending the method of [22] to the Stokes problem. An important feature 
of this method (based on the XFEM approach) is that the Lagrange multiplier is 
identified with the normal trace of the Cauchy stress tensor a(u,p)n at the interface. 
With the aid of stabilization, we have thus a good numerical approximation of this 
quantity, that is crucial in FSI since it gives the force exerted by the viscous fluid 
on the structure. By the way we note that alternative methods based on the work 
of Nitsche [33] (such as [5, 8] in the context of the Poisson problem and [27] in the 
context of the Stokes problems) do not introduce the Lagrange multiplier and thus 
do not necessarily provide a good numerical approximation of this force. 
The outline of the paper is as follows. The continuous problem is set in section 2 
in the weak sense, and the functional spaces are given. We recall the corresponding 
variational formulation with the introduction of a Lagrange multiplier to impose the 
boundary condition in the interface. Next in section 3 the fictitious domain method is 
introduced. In particular the discrete spaces are defined and the discrete variational 
problem is studied without the stabilization technique. This latter - which is an aug- 
mented Lagrangian method - is introduced in section 4, and we show that theoretically 
it enables us to recover the convergence for the multiplier associated with the Dirich- 
let condition (see Lemma 4.2). The convergence analysis for the stabilized method is 
given in section 4.2 and optimal error estimates are proved. Section 5 is devoted to 
numerical tests. Rates of convergence are computed with or without stabilization and 
the behavior of the method is studied for different geometric configurations. More- 
over we compare our method with a classical one which uses a boundary-fitted mesh. 



3 



Technical aspects of the implementation are discussed in section 6. Finally in section 
7 we perform simulations in a simplified unsteady case, what gives a glimpse of the 
future perspectives. The conclusion is given in section 8. 

2. Setting of the problem. 

In a bounded domain of M 2 , denoted by O, we consider a full solid immersed in a 
viscous incompressible fluid. The domain occupied by the solid is denoted by S, and 
we denote by T its boundary. The fluid surrounding the structure occupies the domain 
O \ S — J 7 , where S denotes S U dS (see figure 2.1). 



dO 




Fig. 2.1. Domain for fluid and structure. 



We denote by u and p the velocity field and the pressure of the fluid respectively. 
In this paper, we are interested in the following Stokes problem 

-z/Au + Vp = f in F, (2.1) 

div u = in J 7 , (2.2) 

u = on dO, (2.3) 

u = g on T, (2.4) 

where f G L 2 (J r ), g e H 1/<2 (r). The boundary conditions on V is nonhomogeneous. 
The homogeneous Dirichlct condition we consider on dO has a physical sense, but 
can be replaced by a nonhomogeneous one, without more difficulty. 
With regard to the incompressibility condition, the boundary datum g must obey 

g • ndr = 0. 

r 

We consider this nonhomogeneous condition as a Dirichlet one imposed on T. 
Notice that other boundary conditions are possible on T, such as Neumann conditions, 
as it is done in [22] where mixed boundary conditions are considered. Equation (2.1) 
is the linearized form, in the stationary case, of the underlying incompressible Navier- 
Stokes equations 

du 

— + (u- V)u- i/Au + Wp = f in T . 

The scalar constant v denotes the dynamic viscosity of the fluid. In our presentation, 
for more simplicity, we only consider the stationary case, and the solid is supposed to 
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be fixed. 

The solution of (2.1)-(2.4) can be viewed as the stationary point of the Lagrangian 

L (u,p,\) = v J \D{u)\ 2 dF-J pdiv ud ^-J f • udr - J A • (u - g)dT. 

(2.5) 

Note that we should assume some additional smoothness in (2.5) to make sense, for 
example u e H 2 (J r ), p e H 1 (J r ), A e L 2 (F). The exact solution normally has this 
smoothness provided that f e L 2 (J") and g E H 3 / 2 (F). 

The multiplier A, associated with the Dirichlet condition (2.4), represents the normal 
trace on F of the Cauchy stress tensor. Its expression is given by 

A(u,p) = cr(u,p)n = 2vD(u)n — pn, 

where 

D(u) = \ (Vu + Vu T ) . 

The vector n denotes the outward unit normal vector to dT (see figure 2.1). 

Remark 1. Notice that if we have the incompressibility condition (2.2), then, 
as a multiplier for the Dirichlet condition on F, considering a(u,p)n is equivalent to 
du 

considering v— pn, as it is shown in [20] or [17]. It is mainly due to the equality 

div (Vu + Vu T ) = Au, when div u = 0. 

A finite element method based on the weak formulation derived from (2.5) does 
not guarantee, a priori, the convergence for the quantity cr(u,p)n in L 2 (r). As it has 
been done in [1, 2], our approach consists in considering an augmented Lagrangian in 
adding a quadratic term to the one given in (2.5), as follows 

L(u,p, A) = L (u,p, \)-lJ\\- a(u,p)n| 2 dr. (2.6) 

The goal is to recover the optimal rate of convergence for the multiplier A. The con- 
stant 7 represents a stabilization parameter (see numerical investigations in section 
5.2). It has to be chosen judiciously. 

Let us give the functional spaces we use for the continuous problem (2.1)-(2.4). 
For the velocity u we consider the following spaces 

V = {v g H\T) | v = on dO) , V = Hj(.F), 
V # = {v G V | divv = in J - } , V* = {ve Hj(J") | div v = in J 7 } . 

The pressure p is viewed as a multiplier for the incompressibility condition div u = 0, 
and belongs to L 2 (J r ). It is determined up to a constant that we fix such that p 
belongs to 

Q = LoW = {pe L 2 (^") | ^pdF = oJ. 
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The functional space for the multiplier is chosen as 

w = n-^ 2 {r) = (h 1 / 2 ^))'. 

Remark 2. If we want to impose other boundary conditions, as in [22] for 
instance, the functional spaces V and H 1 / 2 (r) must be adapted, but there is no par- 
ticular difficulty. 

The weak formulation of problem (2.1)-(2.4) is given by: 

Find (u,p, A) G V x Q x W such that 

a(u,v) + 6(v,p) + c(v,A) = £(v) Vv e V, 

b(u,q) = VqeQ, (2.7) 

where 

o(u,v) = 2uJ D(u) : £)(v)dF, (2.8) 

6(u, q) = - f 1 dW udjr ' ( 2 - 9 ) 
c(u,/lx) =~y /x udr, (2.10) 
£(v) = y f ■ vdF, (2.11) 

9{n) = - y^-sdr- ( 2 - 12 ) 

The expression Z?(u) : Z?(v) = trace (D(u)Z)(v) T ) denotes the classical inner product 
for matrices. Let us note that Problem (2.7) is well-posed (see [20] for instance). The 
solution of Problem (2.1)-(2.4) can be viewed as the stationary point of the Lagrangian 
on V x Q x W 

Lq{u,p,\) = v j |D(u)| 2 dJ-- y pdiv ud-^y f ■ u<LF - y A ■ (u - g)dT. 

(2.13) 

3. The fictitious domain method without stabilization. 

3.1. Presentation of the method. 

The fictitious domain for the fluid is considered on the whole domain O. Let us 
introduce three discrete finite element spaces, C H 1 (C) and Q h C Lq(O) on the 
fictitious domain, and W h C L 2 (C). Since O can be a rectangular domain, this 
spaces can be defined on the same structured mesh, that can be chosen uniform (see 
figure 3.1). The construction of the mesh is highly simplified (no particular mesh is 
required). We set 

V h = {v' 1 e C(0) | v\ d0 = 0, vf T G P(T), VT G T h ) , (3.1) 

where P{T) is a finite dimensional space of regular functions such that P(T) 3 Pk(T) 

for some integer k > 1. For more details, see [14] for instance. The mesh parameter 

stands for h = max hr, where is the diameter of T. 
Ter h 
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Fig. 3.1. An example of a mesh on a fictitious domain. The standard degrees of freedom are 
in black, the virtual ones are in red, and the remaining ones are removed. 

Then we define 

V' 1 := Vj^ , Q h := Qf> , W h := Wfc , 

which are natural discretizations of V, L 2 (J r ) and H _1 / 2 (r), respectively. This ap- 
proach is equivalent to XFEM as proposed in [10] or [16] where the standard FEM 
basis functions are multiplied by the Heaviside function 

H(x) - { l -> f ° r X G T 
y > \ 0, for x e O \ T 

and the products are substituted in the variational formulation of the problem. Thus 
the degrees of freedom inside the fluid domain T are used in the same way as in 
the standard FEM, whereas the degrees of freedom in the solid domain S at the 
vertices of the elements cut by the interface (the so called virtual degrees of freedom) 
do not define the field variable at these nodes, but they are necessary to define the 
fields on T and to compute the integrals over T . The remaining degrees of freedom, 
corresponding to the basis functions with support completely outside of the fluid, are 
eliminated (see figure 3.1). We refer to the papers mentioned above for more details. 



Fig. 3.2. Base nodes used for the multiplier space W' 1 . 



An approximation of problem (2.7) is defined as follows 

Find (u h ,p h , \ h ) G V h x Q h x W' 1 such that 

a(u h , v h ) + b{\ h ,p h ) + c(v'\ X h ) = C(v h ) 
b(u h ,q h ) = 

c{M h lt i h ) = g(n h ) 

In matrix notation, the previous formulation corresponds to 




where U, P and A are the degrees of freedom of u h , p h and X h respectively. As it 
is done in [4] or [14] for instance, these matrices ^4„ u , A^ p , A^ x and vectors F°, G° 
are the discretization of (2.8)-(2.12), respectively. Denoting {<£,•}> {xi} and {i/>i} the 
selected basis functions of spaces V , Q h and ~W h respectively, we have 

(A° uu ) .. = 2v f D( Vi ) : D( Vj )&F, (A° up ).. = - f Xi div Vi dF, (A° x ) .. = - f Vi - ^dT, 
(F°). = ^f.^dF, (G°). = -^g-^dr. 

3.2. Convergence analysis. 

Let us define 

V# = |v' 1 g V h | c(v h ,n h ) = 0V/i k e W' 1 } , 
V£ = {v /l G V' 1 | c(v h ,n h ) = c(v h ,g) V/x' 1 € W h } , 
v #,fc = | v ft e yfc | 6( v ' l ,q' 1 ) = Vg h eQ h }, 

V*' h = {v' 1 G V h | &(v'W') = Vq h G Q h ,c(v h ,fi h ) = V/x' 1 G W 1 } . 

The spaces V§, V # - h and Vjf ' h can be viewed as the respective discretizations of the 
spaces V , V# and "Vf. 



Vv' 1 G V' 1 , 

Vq h G Q h , (3.2) 
G W' 1 . 
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Let us assume that the following inf-sup condition is satisfied, for some constant (3 > 
independent of ft: 

HI inf sup „ ,^ V ; g l„ — > p. 

Note that this inf-sup condition concerns only the couple (u,p), and it implies the 
following property 

q h EQ h : b{w h ,q h ) = Vv' 1 G Vq q h = 0. (3.3) 

We shall further assume that the spaces V h , Q h and ~W h are chosen in such a way 
that the following condition is satisfied, for all ft > 

H2 Ji h e W' 1 : c(v\ Ji h ) = Vv h E V h =^JI h = 0. 

Note that this hypothesis is not as strong as an inf-sup condition for the couple ve- 
locity/multiplier. It only demands that the space V h is rich enough with respect to 
the space W h . 

Remark 3. We assume only the inf-sup condition for the couple velocity/pressure, 
not the one for the couple velocity /multiplier . Indeed, the purpose of our work is to 
stabilize the multiplier associated with the Dirichlet condition on T, not the multiplier 
associated with the incompressibility condition. The stabilization of the pressure - on 
the domain T - would be another issue (see page 4^4 of [36] for instance). 

Lemma 3.1. The bilinear form a introduced in (2.8) as 
a : (u,v) ^ 2v [ D(u) : D(v)dJ" 



is uniformly V h -elliptic, that is to say there exists a > independent of ft such that 
for all v h e V h 

a{v h ,v h ) > a\\v h \\ 2 v . 



Proof. Notice that V' 1 C V. Then it is sufficient to prove that the bilinear form 
a is coercive on the space V, that is to say there exists a > such that for all v e V 

a(v, v) > a||v||^. 

By absurd, suppose that for all n € N there exists (v„)„ such that 

«ll-D(v„)ll[L2(^)]4 < ||v„|| v . 

Without loss of generality, we can assume that ||v n ||v = 1- In particular, D(v n ) 
converges to in [L 2 (J r )] 4 . Then, from the Rellich's theorem, we can extract a 
subsequence v m which converges in L 2 (J r ). Using the fact that div v m = 0, the Korn 
inequality (see [14] for instance) enables us to write 

llvm - Vpl&ipr) < C (||v m - VpH^^) + ||D(v ro ) - D{v p )\\ 2 [L2{:F)] ^ , 
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where C denotes a positive constant 1 . This implies that (v m ) TO is a Cauchy se- 
quence in H 1 (J r ). Thus it converges to some which satisfies ||£>(v 00 )|| L 2(-_ F -) = 0. 
The trace theorem implies that we have also v M = on dO. Let us notice that 
v i y ||Z?(v)||[ L 2(jr)]4 is a norm on V. Indeed, if ||D(v 00 )||[L2(_7r)]4 = 0, then is re- 
duced to a rigid displacement, that is to say — I -\- to A x in J~. Then, the condition 
Vqo = on dO leads us to v^, = 0. It belies the fact that ||v m ||v = 1. □ 

Proposition 3.2. Assume that the properties HI and H2 are satisfied. Then 
there exists a unique solution (u h ,p h ,X ) to Problem (3.2). 

Proof. Since Problem (3.2) is of finite dimension, existence of the solution will 
follow from its uniqueness. To prove uniqueness, it is sufficient to consider the case 
f = and g = 0, and to prove that it leads to (u h ,p h ,X h ) = (0,0,0). The last two 
equations in (3.2) show then immediately that \i h £ Vjf'' 1 , so that taking v' 1 = u h in 
the first equation leads to u h = by Lemma 3.1. Taking any test function from Vq 
in the first equation of (3.2) shows now that p h = 0, by condition (3.3) (hypothesis 
HI). And finally the same equation yields X h — by Hypothesis H2. □ 

We recall the following basic result from the theory of saddle point problems [14, 18]. 

Lemma 3.3. Let X and M be Hilbert spaces and A(-, •) : X x X -> R and 
B(-, •) : X x M — > R be bounded bilinear forms such that A is coercive 

A{u,u) > 

and B has the following inf-sup property 

. t B(u,q) 

mi SUp j—r r—r > fj , 

with some a, (3 > 0. Then, for all <fi £ X' and ip e M' , the problem: 

Find u e X and p e M such that 

J A(u,v) + B(v,p) = (<j>,v), \/veX 
\ B(u,q) = (iP,q), VqeM 

has a unique solution which satisfies 

\\u\\x + h\\M<C(U\\ x , + U\\ M ,) 

with a constant C > that depends only on a, ft and on the norms of A and B. 
We can now prove the abstract error estimate for velocity and pressure. 

PROPOSITION 3.4. Assume Hypothesis HI. Let (u,p, A) and (u h ,p h ,X h ) be 
solutions to Problems (2.7) and (3.2) respectively. There exists a constant C > 



1 In the following, the symbol C will denote a generic positive constant which does not depend on 
the mesh size h. It can depend, however, on the geometry of T and T, on the physical parameters, 
on the mesh regularity and on other quantities clear from the context. It can take different values 
at different places. 
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independent of h such that 
Hu-uiv + lb-ZllL^) <C[ inf ||u-v h || v 

1 --h (z\Th 



v'EV 



q h eQ h )j, h ew h J 

(3.4) 



Proof. Take any v h G V^, g' 1 G Q 71 and fi h G W h . Comparing the first lines in 
systems (2.7) and (3.2), we can write 

a(u h -v h ,w h )+b(w h ,p h -q h ) = a(u-v h ,w h )+b(w h ,p-q h )+c(X-fj, h ,w h ) Vw^ G Vft. 

(3.5) 

We have used here the fact that c(A , w' 1 ) = c(fi h ,w h ) = for all w h G Vq . Similarly, 
the second lines in systems (2.7) and (3.2) imply 

b(u h - v h ,s h ) = b(u- v h ,s h ) Vs h eQ h . (3.6) 

Now consider the problem: 

Find x' 1 G V£ and t h G Q h such that 

a(x h ,w h ) + b{w h ,t h ) =a(u- v h ,w h ) + b{w h ,p-q h ) + c(\- fi h ,w h ) Vw h G Vft, 
b{yL h ,s h ) = b(u- v h ,s h ) y S h eQ h . 

Using Lemma 3.3 with A = a, B = b, X = and M = Q h , the solution (x h ,t h ) 
exists and is unique. Moreover, it satisfies 

l|x' l ||v + \\t h \\mr) < C (||u - v' l || v + lb - <z' l || L ^) + II A - ^|| H -V2(r)) ■ 

Comparing the system of equations for (x h ,t h ) with (3.5)-(3.6) and noting that u — 
v' 1 G Vg, we can identify 

x' 1 = u h v h , t h =p h - q h . 
In combination with the triangle inequality, this gives 

||u - u h \\ v + \\p- p h \\ m ^ <C(\\u- v' l || v + \\P q h h H r) + II A - ^Hh-v^d) ■ 

Since G V£, q h G Q h and fi h G W h are arbitrary, this is equivalent to the desired 

result. 

□ 

In summary, the results of this section tell us that, under Hypotheses HI and 
H2, Problem (3.2) has a unique solution which satisfies the a priori estimate (3.4). 
However, we have no estimate for the multiplier X h . 

3.3. The theoretical order of convergence. 

The estimation of the convergence rate proposed for the Poisson problem in [22] can 
be straightforwardly transposed to the Stokes problem. Proposition 3 of [22] ensures 
an order of convergence at least equal to \fh. It can be adapted to our case as follows. 
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Proposition 3.5. Assume Hypotheses HI, H2. Let (u,p,X) be the solution of 
Problem (2.7) for g = 0, such that u E H^p 7 ) n Hq(J") for some e > 0. Assume 
that 

M \\p-q h \\ Q <Ch s , 

q h eQ h 

inf ||A-/i h || w < Ch*, 

li h GW h 

for some S > 1/2. Then 

||u-u h || v + ||p-p fc || La( ^ <CVh. 

Proof. As is shown in [22], Section 3, for any u e H 2+£ (J") n Hq(J") there exists 
a finite element interpolating function w h e Vq such that 

||u-v fc ||v<Cv^. (3.7) 

In fact, v' 1 is constructed as a standard interpolating vector of (1 — %)u where r\h is 
a cut-off function equal to 1 in a vicinity of the boundary T, more precisely in a band 
of width ^r, so that v' 1 vanishes on all the triangles cut by T. This ensures that v' 1 
vanishes on T so that v' 1 <G Vq . Now, the estimate of the present proposition follows 
from (3.4) combined with (3.7) (note that = Vq under our assumptions) and the 
hypotheses on the interpolating functions q h and /j, h . □ 

Let us quote other references that treat of this kind of phenomena, as [18, 37, 38, 
28]. We note, however, that the estimate of the order of convergence in \fh seems too 
pessimistic in view of the numerical tests presented in [22] for the Poisson problem 
(with the possible exception of the lowest order finite elements). In our numerical 
experiments for the Stokes problem, we do not observe the order of convergence as 
slow as \fh. 

4. The fictitious domain method with stabilization. 

4.1. Presentation of the method. 

The main purpose of the stabilization method we introduce consists in recovering the 
convergence on the multiplier A. For that, the idea is to insert in our formulation a 
term which takes into account this requirement. Following the idea used in [1, 2], we 
extend the classical Lagrangian L given in (2.13), as 

L{u,p,\) = vj \D{u)\ 2 &F- j pdiv ud ^-J f-udF-J A (u-g)dr 
-|^|A-a(u,p)n| 2 dr. 

Note that this extended Lagrangian coincides with the previous one on an exact 
solution. The quadratic term so added enables us to take into account an additional 
cost. Minimizing L leads to forcing A to reach the desired value corresponding to 
er(u,p)n. The constant 7 > represents the importance we give to this demand. 
However, notice that this additional term affects the positivity of L. This is the 
reason why we cannot choose 7 too large, and so this approach is not a penalization 
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method. We discuss on this choice of 7 in section 5.2. 
The computations of the first variations leads us to 



5L 
6 



^(v) = 2vj D(u) : D(v)dT - J pdiv vdJ" - J f-v&T-J A • vdr 
+2^7 j (A-<7(u,p)n) • (D(v)n) dT, 
— y qdiv udJ 7 — 7 y g (A — cr(u,p)n) • ndr, 
(/i) = - jT M • (u - g)dr - 7^ (A - a(u,p)n) • /xdT. 
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Thus the stabilized formulation is: 

Find (u,p, A) e V x Q x W such that 
^((u,p,A);v)=£(v) VveV, 



B((u,p,A);g) = Vq e Q, (4.1) 

C((u,p,X);n)=g(n), V/xeW, 



where 



>t((u,p,A);v) = 2u J D(u) : D{v)&F - J pdiv vdJ" - J A • vdr 

-4i/ 2 7 J (L>(u)n) • (L>(v)n) dr + 2t/ 7 ^ p(£>(v)n-n) dr + 2^7 J A • (D(v)n) dr, 
B((u,p,X);q) = - J qdiv udJ"+2^ 7 J q (D(u)n • n) dr - 7 ^ pgdr - 7 ^ gA • ndr, 
C((u,p, A); /lx) = - J V udr + 2^7 y A* ' (-D(u)n)dT -7 J p{n- n)dr J \- y,dT. 
In matrix notation, the previous formulation corresponds to 



^pp Ap\ 

1 uA ^pA ^AA, 



\v A PX I I P I = I 

vAT* A\a / \A/ G, 



where U, P and A are already introduced in section 3.1. As it is done in [4] or [14] 
for instance, these matrices are discretizations of the following bilinear forms 

A uu : (u,v) h- > 2v j D(u) : £>(v)dF - 4i/ 2 7 y (D(u)n) • (D(v)n) dT, 
-4 UP : (v, p) 1— > - y pdiv vdJ" + 2 V1 j p (D(v)n ■ n) dT, 
iuA : (u, A) .— > - y A • vdr + 2^7 y A (D(v)n) dr, 

: (p. 9) 1 — ► -7y P9 dr . 

A>A : (9, A) 1 — >■ -7 y gA • ndr, 
yl A A : ( A, /x) .— > -7 y A • /xdr, 
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and the vectors F and G are the discretization of the following linear forms 

£:v^ J f-vdr, 



g dr. 



Denoting {xi} and {1/^} the selected basis functions of spaces V h , Q ft and W' 1 

respectively, we have 

(A uu ) y - = 2vJ D( Vi ) : DiipJ&F-^iJ WtpJn) ■ {D( Vj n))6T, 
(Ai P )y = - X.div Vi dF + 2i/ 7 jT Xi (I>(Vi)n • n)dr, 

(A»x)« = -7j^Xi(^-- n ) dr ' 
{Axx).. = -1 J MjdT, 

{F).=Jf- Vi dT, {G) l = - J^-^AT. 

4.2. A theoretical analysis of the stabilized method. 

Let us take 7 = 70/1 with some constant 70 > 0. We first observe that the discrete 
problem can be rewritten in the following compact form: 

Find (u h ,p h ,X h ) e V h x Q h x such that 

M{{u h ,p h ,X h );{w h ,q h ,ix h ))=H{w h ,q h ,n h ), V{v h ,q h ,n h ) eV h xQ h x W\ 
where 

M((u,p, A); (v,g,/x)) = 2zv ^ £>(u) : D(\)&F - J (pdivv + gdiv u)dF- J (A • v + n ■ u)dT 
-70/1 y (2i/D(u)n - pn - A) • (2^D(v)n -qn-fi) dT, 



and 



H(v,q,n) = J f • vdr - y /x • gdr. 



In the following, we will need some assumptions for our theoretical analysis: 
Al For all v h e V 1 one has 

/>P(v>|£ 2(r) < C7||v fc |fr. 

A2 For all q' 1 e Q' 1 one has 

h\\q h \\h {r) < C\\q h \\l^ y 
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A3 One has the following inf-sup condition for the velocity-pressure pair of finite 
element spaces 

. b(v h ,q h ) 
^eQ^evJ II^IIl^^IIv^Iv 

with /? > independent of h. 

Assumptions Al and A2 will be discussed in section 5.2 by performing some numer- 
ical tests. 

Note that assumption Al is the same as those introduced in [22] (cf. equations (5.1) 
and (5.5) respectively) in the study of the fictitious domain approach for the Laplace 
equation stabilized a la Barbosa-Hughes. Our assumption A2 is also similar in nature 
to those two, and all these three assumptions can be in fact established if one assumes 
that the intersections of F with the triangles of the mesh are not "too small" (see 
Appendix B of [22] and section 6). Although all these assumptions can be violated 
in practice if a mesh triangle is cut by the boundary Y so that only its tiny portion 
happens to be inside of T . The numerical experiments for the Laplace equation in [22] 
show that such accidents occur rather rarely and their impact on the overall behavior 
of the method is practically negligible. This conclusion can be safely transposed to 
the case of Stokes problem. However, we have now the additional difficulty in the 
form of the inf-sup condition A3. Of course this condition is verified if one chooses 
the classical stable pair of finite element spaces, like for instance the Taylor-Hood ele- 
ments P2/P1 pair for velocity/pressure, and if the boundary P docs not cut the edges 
of the triangles of the mesh. However, in the general case of an arbitrary geometry, 
we have by now no evidence of the fulfillment of the inf-sup condition A3. 

We also need the following result for the L 2 -orthogonal projector from H 1 / 2 (F) 
to W h : 

Lemma 4.1. For allv e H 1/2 (F) one has 

\\P h w - v|| L 2 (r) < CTi 1/2 ||v|| H i/2 ( r), 

where P h denotes the L 2 -orthogonal projector from H 1 / 2 (L) to ~W h . 

Proof. This result is well-known, but we provide for completeness a sketch of the 
proof in the case when discontinuous finite elements are chosen for the space W/j , so 
that W/j contains piecewise constant functions on the mesh on F induced by the 
mesh T h on O (the elements of Tp are the arcs of F obtained by intersecting L with 
the triangles from T' 1 ). The proof in the case of continuous finite elements is similar 
but slightly more technical. 

Let I h be the interpolation operator to the space of piecewise constant functions on 
Tp . For all sufficiently smooth function v on T and for all element tt of the curve F 
obtained by intersection with a triangle T G T h , we set 

I h \\t t = v(x T ), 

where xx is the middle point of tt- We have then I h v e and 

ll^v - v|| L2(r) < ||7 fc v - v|| L2(r) < C7i||v|| H i(r), Vv e H 1 ^), 
by the standard interpolation estimates. Moreover, 



ll^v-v|| L2(r) < ||v|| L2(r) , VveL 2 (r). 
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Interpolating between the last two estimates (see the last chapter of [7]) we get the 
desired result. □ 

We prove in this subsection the following inf-sup result, which is an adaptation of 
Lemma 3 from [22]. 

Lemma 4.2. Under assumptions A1-A3, there exists for 70 small enough a 
mesh-independent constant c > such that 

inf sup M«u\M.^,M)) > c 

(u h ,p h ,x h )ev h xQ h xw h (v', g *,^)ev t xQ'>xW'> |||u /l ,p /l , A III \\\v h , q h , n h \\\ 

where the triple norm is defined by 

\\\u,p, A||| = (\\u\fr + \\p\\l H ^ + h\\D(u)n\\l 2{r) + h\\p\\h ( r) + HM\1 H t) + ^IMIk(r) 



1/2 



and c is a mesh-independent constant. 
Proof. We observe that 

M((u h ,p h , X h ); (u h ,-p h , -X h )) =2i/||u h ||^-7 /i J 4v 2 \D(u h )n\ 2 dT + l0 h J |p fc n - A h | 2 dT 

>^||u' l ||2,+ 7o ^||/n + A' l ||2 2(r) 

where we have used assumption Al and the fact that 70 can be taken sufficiently 
small. More precisely, we can choose 70 such that 4^ 2 7oC < v, where C is the 
constant of assumption Al. The inf-sup condition A3 implies that for all p h € Q h 
there exists v£ G Vg such that 

- jf p"div v^dJ" = \\p h \\h { ^ and ||v£||v < C\\p h \\ h 2 (r) . (4.2) 
Now let us observe that 

M((u h ,p h ,X h );(v h p ,0,0)) = 2u f D(u h ) : D(v h p )dF + \\p h \\h { ^ 

-2v ln h J (2vD(u h )n - p h n - X h ) ■ D(Vp)ndT 
> ll/H^ - „a\\u h \fr - ^\\v% 

-v l0 h a \\2vD(u h )n-p h n-\ h \\i 2{r) ^V(v>|£ 2(r) . 

We have used here the Young inequality which is valid for any a > 0. In particular, 
we can choose a large enough so that we can conclude with the aid of assumptions 
Al and A2 (the constant C here will be independent of a and h, but dependent on 
70 and on the constants in the inequalities Al and A2). We get 

M((u h ,p h ,X h );(v h p ,0,0)) > \\p h \\h { ^ -va\\u h \&- ^-H/llfV) -Cah\\D(u h )n\\l 2{r) 



-^7oMl/n + Al2 2(r) --||/||2 2 
> ^ IbllfV) - Ca\\u h \& v lo ha\\p h n + X h \\l 2(ry 
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Let us now take Ji h = -j i P h \i h where P h is the projector from H 1 / 2 (r) to W h . 
Observe that, in using assumption Al, we have 

M((u h ,p h ,X h ); (0, 0, fi h )) = 1 \\P h u h \\l 2{T) - 7o Jj2vD{u h )rv - p h n - X h ) ■ P h u h dT 

> \\\P h v h \\h { r) 7o (^||£(u>|| L2(r) + Vh\\p h n + \ h \\v (r) ) j= \\P h u h h Hr) 

> ^\\P h " h \\h(r) ~ C\\u h \& Ch\\p h n + \ h \\l H ry 



Combining the above inequalities and taking some small enough numbers k > and 
r] > 0, we can obtain 

M((u h ,p h ,X h ); (u h + nv h p , -p h , X h + V jj, h )) 

> v\\u h \fr + l0 h\\p h n + X h \\l 2{r) + ^\\p h \\h(r) Ca K \\u h \\ 2 v u 7o haK\\p h n + X h \\l 2(r) 
+ ^\\P h u h \\h { r) C V \\n h f v C V h\\p h n + \ h \\l 2{r) 

> %\u h \& + f ll/llfV) + f %>W X h \\l 2{r) + |||PV|| 2 2(r) 

> jWuYv + ^H^ll^r) + ^h\\D(u h )n\\l (T) + 1\\p h \\l (J0 + ^h\\p h \\l (r) + Y h WP hn + A ' l H^(r)- 

In the last line, we have used again assumptions Al and A2 (with the corresponding 
constant C). We now rework the last two terms in order to split p h and \ h . Denoting 

K 

t = — — — , we have 
2Cj 

^h\\P h \\h(r) + | h\\p h n + X h \\l H r) = ? h ((t + l)||/||^ (r) + \\\ h \\h<r) + 2 j/* ■ X h d?j 

> Y h ((* + VWp^Iht) + \\X h \\h<r) (t/2 + l)\\p h \\h (T) ^Vl H^II^Cr)) 

So we finally have 

M((u h ,p h ,X h ); (u h + kv£, X h + m h )) 



> 



c (\\u h \\lr + \\p h \\h (jn + h\\D(u h )n\\l Hr) + h\\p h \\l 2{T) + h\\ X h \\l {r) + ^ll^^ll^(r)) • 



We can now eliminate the projector P h in this estimate by the following calculation, 
which is valid for some (3 > small enough 

|| u *|ft + l -\\P h u h \\l 2{r) > \\u h \& + £ \\P h u h \\l 2{r) = \\u h f v + f (\\u h \\h (T) -\\u h - P h u h \\l 2{r) ) 

> \\u h \fr + £||u fc ||* a(r) Cy?||u"||^ 1/2(r) > \\u h f v + ^\\u h \\h {r) C(3\\u h \\ 2 v 

>l\\n'Tv + l\\u h \\h ( ry 
We have used the result of Lemma 4.1 and the trace inequality. 
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|v£|ft + h\\D(v h p )n\\l Hr) + IWvXht) < C\\v h p \\ 2 Y + C||v£||^ 1/2(r) 

r h\{2 ^ /"yn_/i||2 ^ r<\\\,,h „h \h\ 



In summary, we have obtained that taking 

(v\ q\ » h ) = (U h + KW h p , -p\ X h + V fl h ) 

one has 

M((u h ,p h , \ h ); (v h , q h , fi h )) > c\\\u h ,p h , X h \\\ 2 . (4.3) 
On the other hand, 

11^,^,^111 <M|||u'\/,A' l ||| (4.4) 
with some M > independent of h. Indeed, we have 
|||v h ,g h ,At fc |||< |||u ft ,/,A' i |||+ K |||v^,0,0|||+r,|||0,0,A' l ||| 

< \\\u\p h , \ h \\\ + k (||v£|fr + h\\D(v$)n\\*v (r) + IHWIht)) + r)Vh\\H h L 

Now, by assumption Al and the fact that v£ £ V l so that P h v p l = 0, we have 

Hfv + h\\D(v$)n\\ 2 „ {T) + < C\\v h p \\ 2 v + P^WUvy 

Furthermore, by Lemma 4.1 and by the definition of £ Vq given in (4.2), we have 

<C\WX<C\\p^r)<C\\\n\p\X h \\\. 
We have also 

V^lMii^r) = ^=!|PV l || L2(r) < ±=\\ u h \\v {T) < \\\u\p\X h \\\, 

hence the inequality (4.4). Dividing (4.3) by (4.4) yields 

M{{u h ,p\\ h )-{v\q\y h )) c h h h 
\\\v h ,q h ,H h \\\ ~ M ' P ' 

which is the desired result. □ 

The lemma above, combined with the fact that the bilinear form M is bounded in 
the triple norm on V x Q x W uniformly with respect to h, leads us by a Cea type 
lemma (cf. [14] or Theorem 5.2 in [22]) to the following abstract error estimate 

\\\u-u h , p-p h , \- \ h \\\<C inf |||u- v h ,p-q h ,X- n h \\\. 

(v h ,q h , l l h )£V h xQ h X-W h 

Using the extension theorem for the Sobolev spaces, the standard estimates for the 
nodal (or Clement if necessary) finite element interpolation operators, and the trace 
inequality ||w||L 2 (r) < C ('^IIHIl^t) + ^IIHIl 2 (t)) f° r am/ w € H 1 (T) on any trian- 
gle T &Th (which is valid provided T is sufficiently smooth - see Appendix A of [22] 
for a proof), we obtain the following error estimate 

max(||u - ulv, \\p - p h \\v { r),h\\\ - A h || L2(r) ) < |||u - u h ,p-p h , A - A h ||| 

< C(/i fe "||u|| H ^+i(^) + /i fep+1 ||p!| H fc P +i ( ^) + ^ +1 ||A|| H ^+V2 (r) ), 

where k u , k p and k\ are the degrees of finite elements used for velocity, pressure and 
multiplier A respectively. The proof of this result is rather tedious but can be easily 
reproduced following the ideas of [22] (see, in particular, the proofs of Theorem 5.3 
and Lemma 5.4 there). 



2 (n- 
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5. Numerical experiments. 

For numerical experiments, we consider the square [0, 1] x [0, 1] and choose as T 
the circle whose level-set representation is 

(x - 0.5) 2 + (y - 0.5) 2 = i? 2 , 

with R = 0.21 (see figure 3.1). The exact solutions are chosen equal to 

/ cos(7rx) sm(ny) 



Uex ^'^ ~ ^ - sin(7ra;) cos(Try) 
Pex(x, y) = (y- 1/2) cos(2ttx) + (x - 1/2) sin(2wy). 

The meshes and all the computations have been obtained with the C++ finite 
element library Getfem++ [39]. In the numerical tests, we compare the discrete solu- 
tions with the exact solutions for different meshes (six imbricated uniform meshes). 
We denote U ex , P ex and A ex the discrete forms of functions u ex , p ex and X ex = 
a(u ex ,p ex )n respectively. For practical purposes, the error introduced by the approx- 
imation of the exact vector A ex by A is given by the square root of 

\\A ex - A||L2 (r) = J \cr(U ex , P ex )n - A\ 2 dr. 

This scalar product is developed and using the assembling matrices we compute 

|| Aex - AH^r) = (A uu U ex , U ex ) + 2(A up P ex , U ex ) + 2(A uA A, U ex ) + 

2(A pX A,P ex ) + (A XX A,A), 

where (■, •) denotes the classical Euclidean scalar product in finite dimension. Then, 
the relative error is given by 

\\A ex - A|| L 2 (r) \\A ex - A|| L 2 (r) 

||-^-ea;||L 2 (r) ( ( A UU U ex , U ex) + ( A pp P ex , P ex ) + 2 ( A up P ex , U ex ) ) ^ 

5.1. Numerical experiments for the method without stabilization. 

We present numerical computations of errors when no stabilization arc imposed. We 
consider several choices of the finite element spaces V h , Q h and W h . Four couples 
of spaces are studied (for u/p/A), P1+/P1/P0 (a standard continuous PI element 
for u enriched by a cubic bubble function, standard continuous PI for the pres- 
sure p and discontinuous P0 for the multiplier A element on a triangle), P2/P1/P0, 
for triangular meshes and Q1/Q0/Q0, Q2/Q1/Q0 for quadrangular meshes. The 
elements chosen between velocity and pressure are the ones which ensure the dis- 
crete mesh-independent inf-sup condition HI in the case of uncut functions (ex- 
cept for the Q1/Q0 pair), that is to say the classical case where regular meshes 
are considered. Low degrees are selected to control the memory (CPU time) which 
plays a crucial role in numerical simulations for fluid-structure interactions, spe- 
cially in an unsteady framework. For the multiplier introduced for the interface, 
since the stabilization is not used, a discrete mesh-independent inf-sup condition 
must be satisfied. For instance, the couple of spaces Q1/Q0/Q0 does not satisfy 
this condition. The error curves between the discrete solution and the exact one 
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are given in figure 5.1 for different norms. The rates of convergence are reported. 





Fig. 5.1. Rates of convergence without stabilization for some couples of finite element spaces. 
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The convergence for the fluid velocity is highlighted, whereas the convergence for 
the multiplier seems to not occur, in all cases. We get the convergence for the pressure, 
but not for the test Q1/Q0/Q0 which does anyway not satisfy the inf-sup condition. 
The rates of convergence are better than what we can expect by the theory for u and 
p. The results are not so good for the multiplier. Indeed, without stabilization, the 
the order of magnitude for the relative errors lets us think that the multiplier is not 
well computed. 



5.2. Numerical experiments with stabilization. 

In this part, we consider the method with stabilization terms. Additional terms 
depending on the positive constant 7 are considered in the variational formulation 
(4.1). In the following, we fix 7 = /170, as it is suggested in the proof of Lemma 4.2 
(7 is supposed to be constant, which is natural when uniform meshes are considered). 
The parameter 7 (or 70) has to respond to a compromise between the coercivity of the 
system and the weight of the stabilization term. First, the choice of 7 is discussed. We 
choose the P2/P1/P0 couple of spaces with the space step h = 0.025. To characterize 
a good range of values, we present the condition number (of the whole system) in 
figure 5.2, and the relative errors on the multiplier A for 70 6 [10~ 14 ; 10 4 ] and more 
precisely for 70 € [0.001; 0.200] in figure 5.3. 




Fig. 5.2. The condition number for 70 e [10 14 ;10 4 ]. 
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Fig. 5.3. The relative errors ||A - A h || L 2 (r) for 70 S [10~ 14 ;10 4 ] (left), 70 £ [0.001; 0.200] (right). 

The condition number given for some very small 70 corresponds to the condition 
number of the system when no stabilization is used. For all situations, the condition 
number is degraded when stabilization terms are considered and can explode when 70 
is too large. With regard to the errors on the multiplier A, there is no improvement 
for the relative errors on the multiplier when 70 is too small. When 70 increases, the 
errors on the multiplier becomes interesting even if some peaks can appear (transition 
zone where the coercivity property is very poor). Similar observations (same values 
for 70) are observed on the relative errors for the velocity. 

With regard to the previous experiments, in the following, we choose 70 = 0.05 (so 
7 = 0.05 x h) and we study the numerical convergence analysis of the method when 
stabilization is used. The following numerical experiments have been made in the 
same conditions as the one given in section 3. The results are reported in figure 5.4. 
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Fig. 5.4. Rates of convergence with stabilization for some couples of finite element spaces. 

We notice that we do not observe substantial differences on the rates of conver- 
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gence for the errors on the fluid velocity. As regards to the pressure, a better behavior 
(compared to the first method without stabilization) is observed for the couple of 
spaces Q1/Q0/Q0 that do not satisfy the inf-sup condition. In all cases, the improve- 
ments appear for the multiplier. The method enables to recover the convergence for 
the multiplier. 

5.3. Tests for different geometric configurations. 

In a framework where the solid moves in the fluid domain, we need to perform com- 
putations for different geometric configurations, in order to underline the interest of 
the stabilization method when different types of intersection between the level-set and 
the regular mesh can be achieved. For that, we compute the L 2 (r) relative errors on 
the multiplier A for different positions of the center of the solid, with or without the 
stabilization technique. The perspective is to anticipate the behavior of the method in 
an unsteady case, and these tests enables us to avoid the complexity of a full unsteady 
problem. 

For h = 0.05 and the finite elements triplet P2/P1/P0, we consider the solid as a 
circle, and we make the abscissa of the center of the circle - denoted by xq - vary 
between 0.5 and 0.7 (with a step equal to 0.0005). The variations of the relative 
error (in %) on A are represented in blue (without stabilization) and in red (with 
stabilization). 

10 8 e , , , , , , , , , , 



10 5 - 

<< 
o 




10° I l I l l l l l I l 

0.5 0.52 0.54 0.56 0.58 0.6 0.62 0.64 0.66 0.68 0.7 

*c 

Fig. 5.5. Behavior of the L 2 (r) relative error on A (in serai-log scale), in red with the stabi- 
lization technique (with 70 = 0.05), in blue without. 

In these tests the relevance of our approach using the stabilization technique is 
highlighted when the intersection between the level-set and the mesh varies. Without 
stabilization the errors are huge in many cases (see the curve in blue) , whereas the ro- 
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bustness of the stabilization technique is demonstrated with regards to the constancy 
of the relative errors (see the curve in red). 

5.4. Comparison with a boundary-fitted mesh. 

For three different values of h and by using the elements P2/P1/P0, we compute the 
different relative errors (in %) by using our method (with and without the stabiliza- 
tion technique) and by using a classical code which uses a standard mesh which fits 
closely the boundaries instead of being cut by the boundary of the solid. The results 
are given in Tables 1, 2, 3. 




Fig. 5.6. Different nonuniform boundary-fitted meshes, for which the triangles are not cut. 



h 


L 2 error on u 


H 1 error on u 


L 2 error on p 


L 2 error on A 


0.0358201 
0.0152703 
0.0066282 


0.146643 
0.00371624 
0.00035697 


1.56629 
0.117115 
0.0227257 


3.86771 
0.751358 
0.187311 


9.61603 
3.67841 
1.85277 



Table 1. Errors for a standard uncut mesh. 



h 


L 2 error on u 


H 1 error on u 


L 2 error on p 


L 2 error on A 


0.036418 
0.0150695 
0.00662145 


0.0353448 
0.00274948 
0.00024883 


0.649583 
0.123396 
0.0276422 


2.59781 
0.662703 
0.119263 


6.76061 
13.9277 
1.57377 


Table 2. Errors for a regular cut mesh, without stabilization. 


h 


L 2 error on u 


H 1 error on u 


L 2 error on p 


L 2 error on A 


0.036418 
0.0150695 
0.00662145 


0.03485 
0.00282232 
0.000251731 


0.644208 
0.12423 
0.0275953 


2.46321 
0.556228 
0.104131 


6.61553 
3.71191 
1.52906 



Table 3. Errors for a regular cut mesh, with stabilization (70 = 0.05). 



The results obtained above show that our method enables us to get back the 
precision provided by a classical boundary-fitted mesh. With regards to the errors on 
the multiplier A, notice that by using our method we need to perform the stabilization 
technique in order to recover a good approximation of this variable. 
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5.5. Discussion of assumptions Al and A2. 

In regard to the assumptions Al and A2 considered for the proof of Lemma 4.2, let 
us also study the behavior of the constant C of these assumptions with respect to the 
geometric configuration. In order to verify numerically A2 for instance, we want to 
solve the optimization problem 

h\\q h \\ 2 LHr) h(q h ,q h ) LHr) 
max „ — — = max — -r — rr — . 

One easily shows that the maximum is achieved on the eigenvector of the problem 

%I\xV(r) = A^.xV^) V X h G Q h 

corresponding to the maximal eigenvalue A; = A max . In matrix terms this is rewritten 
as 

^L 2 (r)9i = KA^^q 1 } ^=> hA^j^Aifl^Qi = Kq 1 }, 

where A L 2( r ) and A L 2(jr) are the mass matrices associated with the scalar products 
in L 2 (r) and L 2 (J r ) respectively (see below). Hence the optimal constant in A2 can 
be calculated as X-m^^hA^l^A^^). same thing can be done for Al. Thus we 
consider the two following quantities 

C u (h) — A max (/iA~1^^4 L 2(r)), C p (h) — A max (/i^4L2 1 (jF-)^L 2 (r)), 

where A^rn, Ajji(^), ^L 2 (r)i ^l 2 (.f) denote the matrices respectively defined by 

(^L»(r)) y = I r D(<Pi) ■ D(<Pj)dT, (A m ^).. = Jj.Vtpi : V^cLF + • Vj .<LF, 
( A L 2 (r)) y = Jr Xi ■ Xjdr, (^l 2 (^)) 13 = J^Xi • XjdJ 7 . 

For the particular configuration corresponding to xc = 0.500, let us analyze the 
behavior of max(C u (/i), C p {h)) when the space step h varies. 

300 1 1 1 1 1 1 1 1 , 1 
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Fig. 5.7. Numerical illustration of assumptions Al and A2; max(C u , C p ) in function of h. 
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This graph lets us think that the quantities C u and C p are not constant with 
respect to h (specially when h becomes small), and thus the assumptions Al and A2 
are not satisfied in practice. However, concerning the value of h for which they are 
not satisfied, we get numerically the convergence on the multiplier. At this stage we 
need to consider these assumptions only for proving the theoretical convergence of the 
stabilization technique (see Lemma 4.2). 



6. Some practical remarks on the numerical implementation. 

The numerical implementation of the method for Stokes problem is based on the 
code developed under GetfemH — h Library [39] for Poisson problem. The system 
is solved using the library SuperLU [13]. The advantages of using the GetfemH — h 
library (besides its simplicity of developing finite element codes) is that several specific 
difficulties have been already resolved. Notably, 

- to define basis functions of from traces on T of the basis functions of 
W' 1 . Indeed, their independence is not ensured and numerical manipulations 
must be done in order to eliminate possible redundant functions (and avoid 
to manipulate singular systems), 

- to localize the interface between the fluid and the structure, a level-set func- 
tion which is already implemented (as it is done in [46] for instance), 

- to compute properly the integrals over elements at the interface (during as- 
sembling) external call to Qhull Library [50] is realized (see figure 6.1). 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Fig. 6.1. Local treatment at the interface using Qhull Library. 



As mentioned in the paper [22], it is possible to define a reinforced stability to 
prevent difficulties that can occur when the intersection of the solid and the mesh 
over the whole domain introduce "very small" elements. The technique is based on 
a strategy to select elements which are better to deduce the normal derivative on T. 
A similar approach is given in [35]. This method has been tested for the Dirichlet 
problem in [22], but it is not observed substantial improvements with this enriched 
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stabilization, compared to the results obtained with the stabilization method detailed 
in this paper. However, we expect to take benefits of this second stabilization method 
when the boundary T is led to move through the time, in particular in unsteady 
framework and fluid-structure interactions. 

7. Application to a fluid-structure interaction problem. 

The motivation of our approach lies in the perspective of simulations and control of 
a fluid-solid model for instance. Let us give a simple illustration of that. 

7.1. Coupling with a moving rigid solid. 

In this section, we consider a moving rigid solid which occupies a time-depending 
domain S(t). The displacement of a rigid solid is given by 



where h(t) denotes the coordinates of the center of mass of the solid, and R(i) is 
the rotation which describes the orientation of the solid with respect to its reference 
configuration. In dimension 2, this orientation can be given by a single angle 9(t), 
and we have 



In dimension 2, the angular velocity ui(t) = 9'{t) is a scalar function. The fluid domain 
is given by 0\S(t) = F{t). The state of the corresponding full system is then defined 
by the fluid velocity and pressure, u and p, and the position of the solid given by 
the coordinates of its center of mass h(t) and its angular velocity uj(t). The coupling 
between the fluid and the structure is mainly made at the interface T, through the 
Dirichlet condition 



and through two differential equations which link the position of the solid and the 
forces that the fluid exerts on its boundary, as follows 



Jr(t) 

The vector g denotes the gravity field. Thus, obtaining a good approximation for 
cr(u,p)n is essential for simulating the trajectories of the solid. 

7.2. Illustration: Free fall of a ball. 

The full model described above would necessitate particular attention to the time 
discretization. Indeed, for instance the value of the velocity that we would have to 
consider in the fluid region released by the solid between two time steps has to be 
discussed. Thus, instead of considering the full problem, let us consider a simplified 
approach where the time-dependence aspect is governed only by the position of the 
solid, and not by the time-derivative of the fluid velocity (which requires to tackle the 



%t)=h(t)+R(t)y, ye 5(0), 
S(t) = h(t) + R(i)$(0), 




u(x, t) = W(t) + w(i)(x - h^))- 1 , x e T(t), 




r(t) 



(7.1) 
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difficulty aforementioned). 

A simple illustration consists in simulating in 2D the fall of a rigid ball submitted to 
the gravity force at low Reynolds number. The state of the fluid is then governed 
by the Stokes system we consider in this paper, and the time discretization is only 
about the dynamics of the solid. The radius of the ball is still R = 0.21, and its initial 
position given by the center of the ball C — [xc,yc] — [0.5,0.75]. By symmetry, if 
we assume that the initial velocities are null, then the displacement of the ball is only 
vertical. Thus we impose the Dirichlet condition in the fluid-solid interface as being 
only 

u = h', 

and the function h' = (0, h' 2 ) T satisfies (7.1) which is then reduced to the ID differ- 
ential equation 

Mh' 2 '(t) = -a[h(t)] 2 ti 2 (t) - 9.81M, (7.2) 

where a[h(t)] = / er(u,p)ndr, with T(t) = {h(i) + y | y e T(0)} (the subindcx 2 
Jr(t) 

is used for the second component of the vector), and (u,p) is the solution of 

-i/Au + Vp = in J", 
div u = in J 7 , 
u = on dO, 
u=(0,l) T onT. 

Indeed, the functions u and p are linear we respect to h'. We discretize (7.2) with a 
semi-implicit scheme, as follows 

M 



At 



(h' n 2 +1 h' n 2 ) = -a(h n ) 2 W n 2 +1 9.81M, 

^(hs +i -h5)=hr l . 

For the simulation we choose h = 0.0125 for the space step, still 70 = 0.05 for the 
stabilization parameter, the finite elements triplet P2/P1/P0, At = 10~ 4 for the time 
step, v = 1 and M = 0.02. We represent the amplitude of the velocity at different 
moments in figure 7.1. 
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t = 21 t = 31 £ = 41 



t = 48 t = 51 i = 54 

Fig. 7.1. Simulation of the free fall of a ball in a Stokes flow. 

Note that this simulation cannot be carried out without the stabilization tech- 
nique, because in that case the force that the fluid exerts on the solid is not well- 
computed. Note also that the contact between the ball and the floor would necessitate 
a special treatment that we do not develop here. 

8. Conclusion. 

For Stokes problem which is the corner stone of computations in fluid dynamics, we 
have proposed a fictitious domain method based on extended finite element method. 
Dirichlet boundary conditions at the interface is made using Lagrange multiplier. 
Additional stabilization term is used to ensure an inf-sup condition and to obtain an 
optimal convergence of the normal trace of the Cauchy stress tensor cr(u,p)n. The 
mathematical analysis is presented. We have carried out numerical simulations to 
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compare the new method with the classical finite element approximation based on 
uncut mesh and with the same approach without the introduction of the stabilization 
term. Computations of convergence rates have been performed and have especially 
underlined the interest of the stabilization technique in order to compute a good 
approximation of the normal trace of the Cauchy stress tensor. Besides, this stabi- 
lization technique allows a robust behavior of this quantity when the position of the 
solid changes. 

In a near future, we plan to perform simulations in an unsteady framework, by solving 
the incompressible Navier-Stokes equations in a domain where the solid is moving and 
deforming itself. Our method is particularly interesting in fluid-structure problems 
for which the role of the boundary is central, like for instance when the shape of the 
boundary is the unknown of a control problem. 
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